function TraceComparison_BlueFastVsBlueSlow
%TRACECOMPARISON Summary of this function goes here
%   Detailed explanation goes here

global expPath;

ylim_1 = -200;
ylim_2 = 3600;

%% Define starting dir
if ~isempty(expPath)
    startPath = expPath;
else
    disp('Use SetExpPath to indicate location of colony files')
    return
end
load([startPath '\annotation.mat']);

%% Load SCFP colony
load([startPath '\colony-FOV004-003_SCFP3A.mat']);
colony.exp = annotation.exp;
colony.chip = annotation.chip;

tUnit = colony.exp.TimeInt;
pathCellIndx = GetRedundantPathsII(colony);
list = pathCellIndx(21).lista;
[f, df, time, cellDiv] = FormatPathData(colony, list);
time = tUnit*time;
cellDiv = tUnit*cellDiv;
color = [58 219 255]/255;
colorLine = ([58 219 255]-30)/255;
shade1_h = PlotTrace(time, df, cellDiv, color, colorLine);
str1 = {'Fast (t_{50}=6.4 min), SCFP3A'};
text(240, 3000, str1, 'FontSize', 18)

ylim([ylim_1 ylim_2])
xlabel('\fontsize{17}Time (min)')

set(gca,'box','off')
isholdonque = ishold; 
hold on 
ax = axis; 
plot(ax(2)*[1,1],ax(3:4),'k','linewidth',0.5) 
plot(ax(1:2),ax(4)*[1,1],'k','linewidth',0.5) 
if isholdonque == 0 
hold off 
end

current = get(gca,'DataAspectRatio');
current(2) = 4*current(2);
set(gca, 'DataAspectRatio', current);

%% Load mTurq2 colony
load([startPath '\colony-FOV004-001_mTurq2.mat']);
colony.exp = annotation.exp;
colony.chip = annotation.chip;

tUnit = colony.exp.TimeInt;
pathCellIndx = GetRedundantPathsII(colony);
list = pathCellIndx(23).lista;
[f, df, time, cellDiv] = FormatPathData(colony, list);
time = tUnit*time;
cellDiv = tUnit*cellDiv;
color = [3 146 178]/255;
colorLine = [3 146-30 178-30]/255;
shade2_h = PlotTrace(time, df, cellDiv, color, colorLine);
str1 = {'Slow (t_{50}=33.5 min), mTurquoise2'};
text(205, 3000, str1, 'FontSize', 18)
ylim([ylim_1 ylim_2])
set(gca, 'XTickLabel', []);

current = get(gca,'DataAspectRatio');
current(2) = 4*current(2);
set(gca, 'DataAspectRatio', current);

set(gca,'box','off')
isholdonque = ishold; 
hold on 
ax = axis; 
plot(ax(2)*[1,1],ax(3:4),'k','linewidth',0.5) 
plot(ax(1:2),ax(4)*[1,1],'k','linewidth',0.5) 
if isholdonque == 0 
hold off 
end

%% Show both traces in same axes
fig_h = figure('Position', [300, 100, 900, 600]);
ax = subplot(1,1,1);

figure(shade1_h)
h = get(gcf, 'Children');
newh1 = copyobj(h, fig_h);
posnew = get(newh1,'Position');
posnew(2) = -0.05;
set(newh1, 'Position', posnew)

figure(shade2_h)
h = get(gcf,'Children');
newh2 = copyobj(h, fig_h);
posnew = get(newh2,'Position');
posnew(2) = 0.25;
set(newh2, 'Position', posnew)

figure(fig_h)
str1 = {'\fontsize{17}Fluorescence Production Rate (au)'};
text(-0.1, 0.11, str1, 'Rotation', 90)

axis off
set(gcf, 'Color', 'w')

title_h = title('')

close(figure(1))
close(figure(2))

end

function shade_h = PlotTrace(time, df, cellDiv, color, colorLine)

    threshold  = 400;

    shade_h = figure;
    hold on;
    minmaxT = [time(1) time(end)];
    
    dfBW = df;
    dfBW(df>=threshold)=1;
    indx = df<threshold;
    dfBW(indx) = 0;
        
    stats = regionprops(logical(dfBW), 'BoundingBox');

    nPeak = length(stats);
    plot(time, df, 'linewidth', 1.5, 'color', colorLine);
    for j=1:nPeak
        ti = floor(stats(j).BoundingBox(1));
        tf = ti+ceil(stats(j).BoundingBox(3));
        if j==1, ti = ti+1; end
        if j==nPeak, tf = tf-1; end

        %Calculate intersection with threshold
        m = (df(ti+1) - df(ti))/(time(ti+1) - time(ti));
        b = df(ti) - m*time(ti);
        time_i =  (threshold - b)/m;
        
        m = (df(tf+1) - df(tf))/(time(tf+1) - time(tf));
        b = df(tf) - m*time(tf);
        time_f =  (threshold - b)/m;
        
        
        time_ = [time_i time(ti+1:tf) time_f];
        df_ = [threshold df(ti+1:tf) threshold];
        plot(time_, df_, 'linewidth', 2, 'color', colorLine);
        area(time_, df_, threshold, 'FaceColor', color, 'LineStyle', 'none');

    end
    
    stats = regionprops(logical(~dfBW), 'BoundingBox');
    nPeak = length(stats);
    for j=1:nPeak
       
        ti = stats(j).BoundingBox(1);
        tf = ti+stats(j).BoundingBox(3);

        if tf-ti>=2
            x = [ti+0.5 tf-0.5 tf-0.5 ti+0.5]*5;
            y = [threshold threshold 3600 3600];
            h = patch(x, y, 'r', 'LineStyle', 'none');
            set(h, 'FaceColor', [0.5 0.5 0.5])
            alpha(h, 0.30)
        end
    end
    
    x = [time(1)-100 time(end)+100 time(end)+100 time(1)-100];
    y = [0 0 threshold threshold];
    h = patch(x, y, 'r', 'LineStyle', 'none');
    set(h, 'FaceColor', [0.5 0.5 0.5])
    alpha(h, 0.30)
    
    %Plot lines to indicate cell division
    for j=2:length(cellDiv)
        plot(cellDiv(j), -100, '.k', 'MarkerSize', 13)
    end
    line([-100 2000], [1 1]*0, 'Color', 'k', 'LineStyle','-')
    line([-100 2000], [1 1]*threshold, 'Color', 'k', 'LineStyle','--', 'LineWidth', 2)

    xlim(minmaxT);
    xlim([0 850])
    ylim([-200 2000]);
    box on;
    set(gca, 'FontSize', 13);
    
    hold off;

end

function [f, df, time, cellDiv] = FormatPathData(colony, list)

cell = colony.cell;

time = []; limite = [];
cellDiv = [];   
l = []; dl = [];
ls = []; dls = [];
f = []; df = [];
fs = []; dfs = [];

maxRows = 0;
nImages = 0;
for j=1:length(list)
    
    indx = list(j);
    
    rawLength = max(ceil(cell(indx).bottomLim) - floor(cell(indx).topLim));
    if maxRows<rawLength; maxRows = rawLength; end;
    nImages = nImages + length(cell(indx).frameID);
    
    cellDiv = [cellDiv cell(indx).frameID(1)];
    time = [time cell(indx).frameID];
    f = [f cell(indx).flscnc.f];
    df = [df cell(indx).flscnc.dfFilt];
    l = [l cell(indx).size.lFilt];
    dl = [dl cell(indx).size.dlFilt];
    limite = [limite time(end)];
    
    if j==1
        ls = l;
        dls = dl;
        fs = f;
        dfs = df;
    else
        %Find mother of indx
        indxM = cell(indx).mother;
        %Find sister of indx
        indxS = cell(indxM).daughters(1);
        if indxS == indx, indxS = cell(indxM).daughters(2); end
        
        %Empezar el trace en donde empieza el trace de la hermana
        traceLength = length(cell(indx).flscnc.rawf);
        CurrentTraceLength = length(f);
        indxStartTraceS = CurrentTraceLength - traceLength+ 1;
        
        indxEndTraceS = indxStartTraceS + length(cell(indxS).flscnc.rawf)-1;
        
        ls(indxStartTraceS:indxEndTraceS) = cell(indxS).size.lFilt;
        dls(indxStartTraceS:indxEndTraceS) = cell(indxS).size.dlFilt;
        fs(indxStartTraceS:indxEndTraceS) = cell(indxS).flscnc.f;
        dfs(indxStartTraceS:indxEndTraceS) = cell(indxS).flscnc.df;
    end
end
f = f(1:end-1);
df = df(1:end-1);
time = time(1:end-1);

end
